function [nIters, errRel] = testBlockMG_node3_my(n, TOL, maxIt)
% A driver file to illustrate the usage of block-MG for Michael
%
% INPUT
%  n     : the number of cells, must be of the form 2^m.
%  TOL   : relative error threshold to stop the iteration.
%  maxIt : maximum number of iterations
%
% OUTPUT
%  nIters : the number of iterations performed
%  errRel : the relative error of the solution.
%
% example:
% [nIters, errRel] = testBlockMG_node(1024) should yield the following:
%   iter #1: max-norm of Residual=0.034663; max-norm of error=2.8076e-06
%   iter #2: max-norm of Residual=0.0011727; max-norm of error=4.3362e-08
%   iter #3: max-norm of Residual=4.1648e-05; max-norm of error=1.5591e-09
%   iter #4: max-norm of Residual=1.7096e-06; max-norm of error=6.2258e-11
%   iter #5: max-norm of Residual=7.0835e-08; max-norm of error=2.4709e-12


% a label for using Michaels Matrix, rhs and exact solution
MicOp = 1; %use  = 0 to use the original code here

% This subroutine assembles the Laplacian operator with homogeneous Dirichlet BC.
    function Lap = LaplacianOperator(n, Dim, h)
        % return the 2nd-order Laplacian operator
        %   with homogenuous Dirichlet BC.
        % assuming cell-centered variables and the vector is columnwise
        
        L   = sparse([],[],[],n,n);
        nm1 = ones(n-1,1);
        L   = L + diag(-2*ones(n,1), 0) + diag(nm1, 1) + diag(nm1, -1);
        L = L/h/h;
        
        nP = n^Dim;
        if Dim==1
            Lap = L;
        elseif Dim==2
            In  = sparse(1:n, 1:n, ones(n,1), n, n);
            Lap = kron(L, In) + kron(In, L);
        else
            error('Unsupported value of Dim!');
        end
    end

if nargin<3
    maxIt = 30;
end
if nargin<2
    TOL = 1.e-6;
end

% the label for node-centering.
%  nVoffset==+1 : the boundary nodes are variables;
%  nVoffset==-1 : the boundary nodes are NOT variables.
nVoffset = +1;

addpath('./core','./MicOp');
global mg_matrix;
global mg_param;
mg_param.nu1 = 1;            % number of pre smooths
mg_param.nu2 = 1;            % number of post smooths
mg_param.minlevel = 1;       % level of coarsest grid
mg_param.bx  = 5;            % block size x           
mg_param.by  = 5;            % block size y
mg_param.sx  = 4;            % skip size x
mg_param.sy  = 4;            % skip size y
cycleType = 'f';             % MG cycle type.
                             % the other choice is 'v'
% rectangular grids with node-centered unknowns.
% n is the number of cells.
bds = [-2, 2];
h    = 1/n;
x1d  = linspace(bds(1), bds(2), n+1);
% the following takes the internal nodes as variables.
x1dc = x1d(2:end-1); 
if nVoffset==1
    [X,Y] = meshgrid(x1d, x1d);
elseif nVoffset==-1
    [X,Y] = meshgrid(x1dc, x1dc);
else
    error('Unsupported value of nVoffset!');
end

% some RHS function of the test in page 64 of A Multigrid Tutorial
sinX = sin(X);
sinY = sin(Y);
RHS = -3*sinX.*sinY;
rhs = RHS(:);

% the depth of v-cycle. Again, the number of cells has to be a power of 2.
q=log2(n);
% n is now the number of variables.
n=n+nVoffset;

scaling = 0;
if scaling == 1
    rhs = h*h*rhs; % scale the RHS with h^2 and assign h=1 for the Laplacian operator.
    L2d = LaplacianOperator(n, 2, 1); % assemble the LaplacianOperator
    I = h*h*speye(size(L2d));
else
    %rhs = rhs; % scale the RHS with h^2 and assign h=1 for the Laplacian operator.
    L2d = LaplacianOperator(n, 2, h); % assemble the LaplacianOperator
    I = speye(size(L2d));
end


% change the LaplacianOperator for the fact that
%  we know the values of the unknowns at the boundaries.
L2d = (I - L2d);
change2Identity = 1;
if change2Identity==1 && nVoffset==+1
    % the first column
    L2d(1:n,1:n) = eye(n);
    L2d(1:n,n+1:end) = 0;
    % the upper and lower boundary
    for i=2:n-1
        vid = (i-1)*n+1;
        L2d(vid,:) = 0;
        L2d(vid,vid) = 1;
        vid = i*n;
        L2d(vid,:) = 0;
        L2d(vid,vid) = 1;
    end
    % the last column
    rid = n*n-n+1:n*n;
    L2d(rid,:) = 0;
    L2d(rid,rid) = eye(n);
end

if MicOp %use Michael things instead of things calculated above   
    coeffs = struct('a',1,'b', 1);
    sigma = 1;
    [L2d,rhs,sln] = MicOpWrapper(n,sigma, coeffs);
    
    % sanity check (assertion will fall if the condition doesn't sattisfied)    
    assert( all(L2d\rhs - sln < 0.01))
else
     % the exact solution we expect.
    sln = L2d\rhs;
end

% set the gobal matrices and initialize the MG 
mg_matrix = mg_matrix_init(L2d, q, 'p1', nVoffset);

% use zero as the initial guess.
uS = zeros(n*n,1); 
% the initial residual
res0 = norm(rhs,inf);
% the initial error
err0 = norm(sln,inf);
for jj=1:maxIt
    if cycleType=='v'
%        uS = vcycle_block(uS, rhs, q);
        uS = vcycle(uS, rhs, q);
    elseif cycleType=='f'
%        uS = fcycle_block(uS, rhs, q);
        uS = fcycle(uS, rhs, q);
    else
        error('unknown cycle type!');
    end

    tmp = rhs-L2d*uS;
    res = norm(tmp,inf);
    resRel = res/res0;
    errRel = norm(uS-sln, inf)/err0;
    
    info = ['  iter #', int2str(jj),'  Grid Size=', int2str(n),  ...
        ': max-norm of Residual=', ...
        num2str(resRel,'%5.2e'), '; max-norm of error=', num2str(errRel,'%5.2e')];
    disp(info);
%    if (errRel<TOL)
    if (resRel<TOL)
        break
    end
end
nIters = jj;

end


%this function create an operator to the problem 
% a*u_xx+b*u_xx - sigma*u = f 
% on (-2:2)x(-2,2) domain
% with dirichlet boundary condition boundary condition
% u = g at boundary 
% The boundary conditions are imposed inside the operator.
%
% Input:
% N      - the size of grid
% sigma  - the coefficient in the equation above
% coeffs - a structure with constant coefficients a and b in the equation above
%
% Output
% A   - the discrete operator that approximate the problem above
% rhs - a vector rhs of the problem
% sln - a vector exact solution to the problem approimated by linear system Ax=rhs
%
function [A,rhs,sln] = MicOpWrapper(N,sigma, coeffs)

    Grid  = Tools.Grid.CartesianGrid( -2,2,N,-2,2,N);                  
    
    LapOp =  LaplacianOp(Grid,sigma, coeffs,0);

    Exact =  sin(coeffs.a.*Grid.X).*sin(coeffs.b.*Grid.Y);
    
    RHS = -coeffs.a*Exact -coeffs.b*Exact - sigma.*Exact;   
    RHS = LapOp.AdjustRhs(RHS,Exact);

    rhs = RHS(:);   
    sln = Exact(:);    
    A   = LapOp.A;       
end
